Atelier R Cartographie
L’écosystème spatial sur R
Introduction à sf
Site web de
sf: Simple Features for Rsf pour Simple Features
Sortie en octobre 2016
A pour but de rassembler les fonctionnalités d’anciens packages (
sp,rgeosandrgdal) en un seulFacilite la manipulation de données spatiales, avec des objets simples.
Tidy data: compatible avec la syntaxe pipe
%>%et les opérateurs dutidyverse.Principal auteur et mainteneur : Edzer Pebesma (également auteur du package
sp)
la structure des objets sf :
Importer / exporter des données
library(sf)## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
mtq <- read_sf("data/mtq/martinique.shp")
mtq <- st_read("data/mtq/martinique.shp")## Reading layer `martinique' from data source `/home/comeetie/Projets/quantilille/lecture/data/mtq/martinique.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34 features and 23 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
## Projected CRS: WGS 84 / UTM zone 20N
write_sf(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)
st_write(mtq,"data/mtq/martinique.gpkg",delete_layer = TRUE)## Deleting layer `martinique' using driver `GPKG'
## Writing layer `martinique' to data source `data/mtq/martinique.gpkg' using driver `GPKG'
## Writing 34 features with 23 fields and geometry type Polygon.
Système de coordonées
Les projections/système de coordonées sont répertoriées grâce à un code le code epsg :
- lat/long : 4326 https://epsg.io/4326
- Lambert 93 : 2154 https://epsg.io/2154
- Pseudo-Mercator : 3857 https://epsg.io/3857
- Lambert azimuthal equal area : 3035 https://epsg.io/3035
Projection
Obtenir la projection en utilisant st_crs() (code epsg) et la modifier en utilisant st_transform().
st_crs(mtq)## Coordinate Reference System:
## User input: WGS 84 / UTM zone 20N
## wkt:
## PROJCRS["WGS 84 / UTM zone 20N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 20N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-63,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## ID["EPSG",32620]]
mtq_4326 <- mtq %>% st_transform(4326)Afficher les données
Affichage par défaut :
plot(mtq)## Warning: plotting the first 10 out of 23 attributes; use max.plot = 23 to plot
## all
En ne gardant que la géométrie :
plot(st_geometry(mtq))Extraire les centroïdes
mtq_c <- st_centroid(mtq)## Warning in st_centroid.sf(mtq): st_centroid assumes attributes are constant over
## geometries of x
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)Matrice de distance
mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]## Units: [m]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000 35297.56 3091.501 12131.617 17136.310
## [2,] 35297.557 0.00 38332.602 25518.913 18605.249
## [3,] 3091.501 38332.60 0.000 15094.702 20226.198
## [4,] 12131.617 25518.91 15094.702 0.000 7177.011
## [5,] 17136.310 18605.25 20226.198 7177.011 0.000
Agrégation de polygones
Union simple :
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")A partir d’une variable de regroupement :
library(dplyr)##
## Attachement du package : 'dplyr'
## Les objets suivants sont masqués depuis 'package:stats':
##
## filter, lag
## Les objets suivants sont masqués depuis 'package:base':
##
## intersect, setdiff, setequal, union
mtq_u2 <- mtq %>%
group_by(STATUT) %>%
summarize(P13_POP = sum(P13_POP))
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u2), add=T, lwd=2, border = "red", col=NA)Zone tampon
mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")Intersection de polygones
# create a polygon
m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586),
c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)st_intersection() extrait la partie de mtq qui s’intersecte avec le polygone créé.
mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)Compter les points dans des polygones
st_sample() crée des points aléatoires sur la carte de Paris.
pts <- st_sample(x = mtq, size = 50)
plot(st_geometry(mtq))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)st_interects() crée une liste de points dans chaque polygone.
inter <- st_intersects(mtq, pts)Il ne reste plus qu’a compter.
mtq$nbpts <- sapply(X = inter, FUN = length)
plot(st_geometry(mtq))
# display munucipalities that intersect at least 2 point
plot(st_geometry(mtq[mtq$nbpts>2,]), col = "grey", add=TRUE)
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)Autre solution, faire une jointure spatiale et aggréger !
mtq_counts <- mtq %>% st_join(st_as_sf(pts)) %>% count(INSEE_COM)
plot(mtq_counts %>% select(n))
plot(pts, pch = 20, col = "red", add=TRUE, cex = 1)Voronoi Polygons
google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)
mtq_v <- st_collection_extract(st_voronoi(x = st_union(mtq_c)))
mtq_v <- st_intersection(mtq_v, st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c)
plot(st_geometry(mtq_v), col='lightblue')Autre traitements
- st_area(x)
- st_length(x)
- st_disjoint(x, y, sparse = FALSE)
- st_touches(x, y, sparse = FALSE)
- st_crosses(s, s, sparse = FALSE)
- st_within(x, y, sparse = FALSE)
- st_contains(x, y, sparse = FALSE)
- st_overlaps(x, y, sparse = FALSE)
- st_equals(x, y, sparse = FALSE)
- st_covers(x, y, sparse = FALSE)
- st_covered_by(x, y, sparse = FALSE)
- st_covered_by(y, y, sparse = FALSE)
- st_equals_exact(x, y,0.001, sparse = FALSE)
- …
Conversion
- st_cast
- st_collection_extract
- st_sf
- st_as_sf
- st_as_sfc
Autres packages
CRAN task views permet d’avoir des informations sur les packages du CRAN pertinents pour des tâches reliées à certains sujets.
CRAN Task View: Analysis of Spatial Data:
- Classes for spatial data
- Handling spatial data
- Reading and writing spatial data
- Visualisation
- Point pattern analysis
- Geostatistics
- Disease mapping and areal data analysis
- Spatial regression
- Ecological analysis
Préparer / récupérer des données
Premier exemple les données sont stockées dans un fichier shape file. Tout devrait bien ce passer:
library(sf)
library(dplyr)
# Import de la couche géographique (iris de Paris)
iris.75 <- st_read("../data/iris_paris.shp", stringsAsFactors = F) ## Reading layer `iris_paris' from data source `/home/comeetie/Projets/quantilille/data/iris_paris.shp' using driver `ESRI Shapefile'
## Simple feature collection with 992 features and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
## Projected CRS: RGF93 / Lambert-93
Regardons la projection utilisée:
st_crs(iris.75)## Coordinate Reference System:
## User input: RGF93 / Lambert-93
## wkt:
## PROJCRS["RGF93 / Lambert-93",
## BASEGEOGCRS["RGF93",
## DATUM["Reseau Geodesique Francais 1993",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4171]],
## CONVERSION["Lambert-93",
## METHOD["Lambert Conic Conformal (2SP)",
## ID["EPSG",9802]],
## PARAMETER["Latitude of false origin",46.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8821]],
## PARAMETER["Longitude of false origin",3,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8822]],
## PARAMETER["Latitude of 1st standard parallel",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8823]],
## PARAMETER["Latitude of 2nd standard parallel",44,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8824]],
## PARAMETER["Easting at false origin",700000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8826]],
## PARAMETER["Northing at false origin",6600000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8827]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["France - onshore and offshore, mainland and Corsica."],
## BBOX[41.15,-9.86,51.56,10.38]],
## ID["EPSG",2154]]
Un autre cas de figure fréquent, des données ponctuelles stockées dans un csv avec deux colonnes ( latitude et longitude en wgs 84). Dans ce cas, on importe le csv puis on convertit la data.frame en data.frame spatiale (sf) avec la fonction st_as_sf. Il suffit de spécifier le nom des colones contenant les coordonnées et le CRS, généralement on utilisera le code epsg de wgs 84 (4326).
# Import du dataset
accidents.2019.paris <- readRDS("../data/accidents2019_paris.RDS")
# Transformation en objet sf
accidents.2019.paris <- st_as_sf(accidents.2019.paris,
coords = c("long", "lat"),
crs = 4326, agr = "constant") %>%
st_transform(2154)
plot(st_geometry(accidents.2019.paris))Regardons rapidement ces deux jeux de données et comptons le nombre de blessés graves par iris.
inter <- st_intersects(iris.75, accidents.2019.paris)
inter_blessgravtues <- st_intersects(iris.75, accidents.2019.paris
%>% filter(grav%in%c(2,3)))
iris.75$nbacc <- sapply(X = inter, FUN = length)
iris.75$nbacc_blessgravtues <- sapply(X = inter_blessgravtues, FUN = length)
#Il manque 24 accidents
nrow(accidents.2019.paris)-sum(iris.75$nbacc)## [1] 24
plot(st_geometry(iris.75))
plot(accidents.2019.paris %>% filter(grav%in%c(1,4)) ,
pch = 20, col = "darkgreen", add=TRUE, cex = 0.5)## Warning in plot.sf(accidents.2019.paris %>% filter(grav %in% c(1, 4)), pch =
## 20, : ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav%in%c(3)) ,
pch = 20, col = "orange", add=TRUE, cex = 0.5)## Warning in plot.sf(accidents.2019.paris %>% filter(grav %in% c(3)), pch = 20, :
## ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav==2),
pch = 20, col = "red", add=TRUE, cex = 1)## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 2), pch = 20, :
## ignoring all but the first attribute
Utiliser osmdata
osmdata permet d’extraire des éléments de la base de données gratuite et open-source OpenStreetMap. Cela peut nous servir par exemple pour récupérer des élements d’habillage fleuve / routes ou d’autres informations. La requête suit la nomenclature osm sur base de clés / valeures. Vous pouvez utiliser tagingo pour explorer l’ensemble des clés et valeures utilisées par la communauté OSM.
library(osmdata)
# Récupérer les routes principales grâce à osm
bb <- iris.75 %>% st_transform(4326) %>% st_bbox()
q <- opq(bbox = bb,timeout = 180)
qm <- add_osm_feature (q, key = 'highway',value = 'motorway', value_exact = FALSE)
qt <- add_osm_feature (q, key = 'highway',value = 'trunk', value_exact = FALSE)
qp <- add_osm_feature (q, key = 'highway',value = 'primary', value_exact = FALSE)
motorway<- osmdata_sf(qm)
trunk <- osmdata_sf(qt)
primary <- osmdata_sf(qp)
roads <- c(primary,trunk,motorway)$osm_lines %>% st_transform(st_crs(iris.75))
roads.geom = st_geometry(roads)
# Récupérer le shape de la seine
qr <- q %>%
add_osm_feature (key = 'waterway') %>%
add_osm_feature(key = "name:fr", value = "La Seine")
river <- osmdata_sf(qr)
river.geom <- c(st_geometry(river$osm_lines),st_geometry(river$osm_multilines)) %>% st_transform(st_crs(iris.75))
# Export road and river layers to shapefile
st_write(roads.geom, dsn = "data/osmdata/road.shp")
st_write(river.geom, dsn = "data/osmdata/river.shp")Utilisons ces données pour habiller un peu notre carte:
# bbox est utilisé pour centrer sur Paris
bb <- st_bbox(iris.75)
par(mar=c(0.2,0.2,1.4,0.2), bg="ivory")
plot(st_geometry(iris.75), col = "ivory", border="ivory3",
xlim = bb[c(1,3)], ylim = bb[c(2,4)])
plot(st_geometry(roads.geom),col="#666666",lwd = 1.2,add=TRUE)
plot(st_geometry(river.geom),col="#87cdde",lwd = 3,add=TRUE)
plot(accidents.2019.paris %>% filter(grav==3) , pch = 20, col = "orange", add=TRUE, cex = 1)## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 3), pch = 20, :
## ignoring all but the first attribute
plot(accidents.2019.paris %>% filter(grav==2) , pch = 20, col = "red", add=TRUE, cex = 1)## Warning in plot.sf(accidents.2019.paris %>% filter(grav == 2), pch = 20, :
## ignoring all but the first attribute
Géocodage
Géocoder c’est passé d’une adresse à une position géographique. En france, la Base d’adresse national permet de faire ce travail efficacement, en R le package banR permet d’intéroger l’API de la BAN simplement ce package doit être installé via github et permet ensuite de géocoder une colonne d’adresse en batch en spécifiant la colone contenant les adresse et éventuellement une colone contenant un code insee pour faciliter et préciser la requête. Pour des adresses internationales il est possible d’utiliser tidygeocoder qui peut intéroger différente API gratuite ou payante et qui fonctionne de manière assez similaire.
# avec banR
# remote
library(banR)
geo_banR = accidents.2019.paris %>%
filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>%
geocode_tbl(adresse = voie,code_insee = com) %>% select(latitude,longitude) %>%
st_as_sf(coords = c("longitude", "latitude"),crs = 4326, agr = "constant") %>%
st_transform(2154)## Writing tempfile to.../tmp/RtmpIrPGvn/file239a2d1957f8.csv
## Warning: The `path` argument of `write_csv()` is deprecated as of readr 1.4.0.
## Please use the `file` argument instead.
## If file is larger than 8 MB, it must be splitted
## Size is : 666 bytes
## SuccessOKSuccess: (200) OK
## New names:
## * geometry -> geometry...10
## * geometry -> geometry...13
library(tidygeocoder)
geo_tidygeocoder = accidents.2019.paris %>%
filter(catv %in% c("VAE","EDP à moteur")) %>% slice(1:10) %>%
mutate(addr = paste(voie,", Paris, France")) %>%
geocode(addr,method="osm") %>% select(lat,long) %>%
st_as_sf(coords = c("long", "lat"),crs = 4326, agr = "constant") %>%
st_transform(2154)
st_distance(geo_banR,geo_tidygeocoder,by_element = TRUE)## Units: [m]
## [1] 950.52141 382.66079 122.00239 181.62935 47.91376 642.02833 41.90755
## [8] 305.16959 664.76010 664.76010
Faire des cartes interactive avec R
De nombreuses solutions existent pour faire des cartes interactives avec R :
mapview,leafletetmapdeckpermettent de faire des cartes interactives.
Par simplicité, nous nous concentrons ici sur ggplot2 pour la partie statique et mapview pour la partie interactive.
Petite introduction de sémiologie graphique
Cartes interactives mapview
Les cartes interactives ne sont pas forcément très pertinentes pour représenter des informations géostatistiques.
En revanche, elles sont utiles pour explorer les bases de données. Voyons un exemple avec mapview concernant les accidents mortels à Paris en 2019.
#remotes::install_github("r-spatial/mapview")
library(mapview)
mapviewOptions(fgb=FALSE)
# construire une carte avec certaines options pour les cercles
# avec mapview la taille des cercles reste constante quel que soit le zoom.
# grav = 2 : individus tués
mapview(accidents.2019.paris %>%
filter(grav==2))On customise un peu…
mapview(accidents.2019.paris %>%
filter(grav==2),
map.types = "Stamen.TonerLite", legend=FALSE,
cex=5, col.regions="#217844", lwd=0, alpha=0.9)On customise encore un peu plus…
mapview(accidents.2019.paris %>%
filter(grav==2) %>%
mutate(age=2019-an_nais),
map.types = "Stamen.TonerLite", legend=FALSE,
cex="age", zcol="sexe", lwd=0, alpha=0.9)Faire des cartes statiques
ggplot2est un package très utilisé pour faire tous types de graphiques, et a été adapté spécifiquement aux cartes (geom_sf).- Le package
tmapcontient des fonctionnalités avancées basées sur la logique deggplot2 map_sf(anciennementcartography) s’appuie sur un langage dit “base R” et permet de faire des représentations cartographiques, basiques comme avancées.
ggplot2
grammar of graphics
- “The Grammar of Graphics” (Wilkinson, Annand and Grossman, 2005)
- grammaire → même type de description pour des graphique différents
Composants de la grammaires :
ex : f(data) → x position, y position, size, shape, color
ex : points, lines, bars, texts
ex : f([0, 100]) → [0, 5] px
ex : segmentation des données suivant un ou plusieurs facteurs
ex : moyenne, comptage, régression
Création d’un graphique :
! Données toujours sous forme de data.frame bien formatées
Exemple diagramme en barre :
library(ggplot2)
ggplot(accidents.2019.paris)+geom_bar(aes(x=catv,group=sexe,fill=sexe))Qui mérite quelques ajustements :
catv_ol = accidents.2019.paris %>% st_drop_geometry %>% count(catv) %>% arrange(n) %>% pull(catv)
gg = accidents.2019.paris %>% mutate(catv_o = factor(catv,levels=catv_ol)) %>% filter(catv_o %in% tail(catv_ol,10))
ggplot()+geom_bar(data = gg,aes(y=catv_o,group=sexe,fill=sexe))+
scale_fill_brewer("Sexe",palette="Set1")+
theme_bw()+
labs(title="Nombre d'accidentés par type de véhicule et sexe",
subtitle="à Paris en 2019, pour les hommes et les femmes ",caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")Quelques chose de plus exotique :
gg = accidents.2019.paris %>%
st_drop_geometry %>%
filter(catv %in% tail(catv_ol,9)) %>%
count(catv,lum,sexe) %>%
add_count(catv,wt=n,name="tot") %>%
mutate(prop = n/tot)
ggplot(gg)+geom_point(aes(y=lum,x=sexe,color=prop,size=prop))+
facet_wrap(~catv)+
scale_color_distiller(palette="Reds",direction=1)+
labs(title="Part d'accidentés par type de véhicule et sexe",
subtitle="à Paris en 2019. ",caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")Ou bien encore :
library(tidyr)
catv_sel = c("Bicyclette","VL seul","VAE","VU seul","EDP à moteur","Scooter < 50 cm3")
gg = accidents.2019.paris %>%
st_drop_geometry %>%
filter(catv %in% catv_sel) %>%
count(catv,sexe) %>% pivot_wider(names_from = "sexe",values_from="n")
ggplot(gg) + geom_segment(aes(x="Homme",y=Masculin,xend="Femme",yend=Féminin,color=catv)) ggplot(gg)+
geom_segment(aes(x="Homme",y=Masculin,xend="Femme",yend=Féminin,color=catv))+
geom_text(data=gg %>% filter(catv!="Scooter < 50 cm3"),aes(x="Homme",y=Masculin,label=catv,color=catv),hjust="left")+
geom_text(data=gg %>% filter(catv!="VU seul"),aes(x="Femme",y=Féminin,label=catv,color=catv),hjust="right")+
scale_color_discrete(guide="none")+theme_bw()+
labs(title = "Nombre d'accidentés suivant le sexe et le type de véhicule",
subtitle="à Paris en 2019", caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")Intégrer des données spatiales avec geom_sf
Cartes avec ronds proportionnels
ggplot() +
geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
geom_sf(data = river.geom, colour = "azure",size=2) +
geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
geom_sf(data = iris.75 %>% st_centroid(),
aes(size= nbacc), colour="#E84923CC", show.legend = 'point') +
scale_size(name = "Nombre d'accidents",
breaks = c(1,10,100,200),
range = c(0,5)) +
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(iris.75)[c(1,3)],
ylim = st_bbox(iris.75)[c(2,4)]) +
theme_minimal() +
theme(panel.background = element_rect(fill = "ivory",color=NA),
plot.background = element_rect(fill = "ivory",color=NA)) +
labs(title = "Nombre d'accidents de la route à Paris par iris",
caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021",x="",y="")## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x
Cartes choroplethes
library(RColorBrewer)
acc = iris.75 %>%
st_join(accidents.2019.paris) %>%
group_by(INSEE_COM,do_union=TRUE) %>%
summarize(nb_acc=n(),
nb_vl=sum(if_else(catv=="VL seul",1,0),na.rm = TRUE),
nb_edp=sum(if_else(catv=="EDP à moteur",1,0),na.rm = TRUE),
nb_velo=sum(if_else(catv=="Bicyclette",1,0),na.rm = TRUE))## `summarise()` has grouped output by 'INSEE_COM'. You can override using the `.groups` argument.
acc## Simple feature collection with 20 features and 6 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 643075.6 ymin: 6857477 xmax: 661086.2 ymax: 6867081
## Projected CRS: RGF93 / Lambert-93
## # A tibble: 20 x 7
## # Groups: INSEE_COM [20]
## INSEE_COM do_union nb_acc nb_vl nb_edp nb_velo geometry
## <chr> <lgl> <int> <dbl> <dbl> <dbl> <POLYGON [m]>
## 1 75101 TRUE 285 117 8 27 ((651932 6861789, 651927.1 68…
## 2 75102 TRUE 173 63 10 19 ((650718.1 6863519, 650706.1 …
## 3 75103 TRUE 202 63 10 21 ((652440 6862594, 652433.4 68…
## 4 75104 TRUE 277 80 25 41 ((651966.6 6861729, 651911.6 …
## 5 75105 TRUE 234 75 14 28 ((651960.7 6859931, 651886.2 …
## 6 75106 TRUE 273 112 16 30 ((650807.6 6860419, 650795.2 …
## 7 75107 TRUE 444 202 13 37 ((650112.5 6861125, 650080 68…
## 8 75108 TRUE 721 333 17 32 ((650208.6 6862824, 650181.1 …
## 9 75109 TRUE 320 118 11 41 ((650787.8 6863551, 650780.1 …
## 10 75110 TRUE 442 162 11 57 ((652193.8 6863596, 652167.8 …
## 11 75111 TRUE 509 171 20 69 ((653355.9 6863227, 653352.2 …
## 12 75112 TRUE 1000 441 13 77 ((655285.2 6858619, 655221.2 …
## 13 75113 TRUE 798 350 15 62 ((651662.8 6858694, 651667.5 …
## 14 75114 TRUE 704 266 19 49 ((650148.5 6858070, 650116.1 …
## 15 75115 TRUE 858 364 14 78 ((645885.5 6859571, 645942.5 …
## 16 75116 TRUE 1429 766 16 67 ((645783.9 6859590, 645571.2 …
## 17 75117 TRUE 903 414 12 47 ((648839.5 6864557, 648790.3 …
## 18 75118 TRUE 718 318 17 46 ((651443.5 6864903, 651436.5 …
## 19 75119 TRUE 855 386 16 65 ((653702.1 6865070, 653697.5 …
## 20 75120 TRUE 826 358 15 55 ((654559 6863398, 654461.5 68…
bks <- round(10*quantile(100*acc$nb_velo/acc$nb_acc, na.rm=TRUE, probs=seq(0,1,0.2)))/10
bks## 0% 20% 40% 60% 80% 100%
## 4.4 6.6 7.7 9.8 12.1 14.8
acc = acc %>% mutate(txaccvelo= 100*nb_velo/nb_acc,txaccvelo_cat = cut(txaccvelo,bks))
ggplot() +
geom_sf(data = iris.75,colour = "ivory3",fill = "ivory") +
geom_sf(data = acc, aes(fill = txaccvelo_cat)) +
geom_sf(data = river.geom, colour = "#87cdde",size=2) +
geom_sf(data = roads.geom, colour = "#666666",size=0.5) +
scale_fill_brewer(name = "Part (En %)",palette = "Reds", na.value = "grey80") +
coord_sf(crs = 2154, datum = NA,
xlim = st_bbox(iris.75)[c(1,3)],
ylim = st_bbox(iris.75)[c(2,4)]) +
theme_minimal() +
theme(panel.background = element_rect(fill = "ivory",color=NA),
plot.background = element_rect(fill = "ivory",color=NA)) +
labs(title = "Part des Accidentés à vélos",
subtitle = "par arrondissement à Paris en 2019",
caption = "fichier BAAC 2019, ONISR\nantuki & comeetie, 2021", x="",y="")Crédits et reproductibilité
Présentation faite grâce au package rmdformats.
Elle s’inspire très fortement, ainsi que son tutoriel, d’une précédente formation donnée par les mêmes auteurs avec Timothée Giraud.
Partage de la configuration de R et des packages utilisés :
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.2 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 tidyr_1.1.3 ggplot2_3.3.3 mapview_2.9.9
## [5] dplyr_1.0.6 sf_0.9-8
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 svglite_2.0.0 lattice_0.20-44
## [4] png_0.1-7 class_7.3-19 leaflet.providers_1.9.0
## [7] assertthat_0.2.1 digest_0.6.27 utf8_1.2.1
## [10] R6_2.5.0 leafpop_0.1.0 stats4_4.1.0
## [13] evaluate_0.14 e1071_1.7-7 highr_0.9
## [16] pillar_1.6.1 rlang_0.4.11 uuid_0.1-4
## [19] rstudioapi_0.13 raster_3.4-10 rmarkdown_2.8
## [22] labeling_0.4.2 webshot_0.5.2 stringr_1.4.0
## [25] htmlwidgets_1.5.3 munsell_0.5.0 proxy_0.4-25
## [28] compiler_4.1.0 xfun_0.23 pkgconfig_2.0.3
## [31] systemfonts_1.0.2 base64enc_0.1-3 htmltools_0.5.1.1
## [34] tidyselect_1.1.1 tibble_3.1.2 bookdown_0.22
## [37] codetools_0.2-18 fansi_0.5.0 crayon_1.4.1
## [40] withr_2.4.2 grid_4.1.0 jsonlite_1.7.2
## [43] satellite_1.0.2 gtable_0.3.0 lifecycle_1.0.0
## [46] DBI_1.1.1 magrittr_2.0.1 units_0.7-1
## [49] scales_1.1.1 rmdformats_1.0.2 KernSmooth_2.23-20
## [52] cli_2.5.0 stringi_1.6.2 farver_2.1.0
## [55] leaflet_2.0.4.1 sp_1.4-5 ellipsis_0.3.2
## [58] brew_1.0-6 generics_0.1.0 vctrs_0.3.8
## [61] tools_4.1.0 leafem_0.1.6 glue_1.4.2
## [64] purrr_0.3.4 crosstalk_1.1.1 yaml_2.2.1
## [67] colorspace_2.0-1 classInt_0.4-3 knitr_1.33